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In this paper an SPH version of the alpha turbulence model devised by Holm 
and his colleagues is formulated for compressible flow with a resolution that varies 
in space and time. The alpha model involves two velocity fields. One velocity field 
is obtained from the momentum equation, the other by averaging this velocity field 
as in the version of SPH called XSPH. The particles (fluid elements) are moved 
with the averaged velocity. In analogy to the continuum alpha model we obtain 
a particle Lagrangian from which the SPH alpha equations can be derived. The 
system satisfies a discrete Kelvin circulation theorem identical to that obtained 
with no velocity averaging. In addition, the energy, linear and angular momentum 
are conserved. We show that the continuum equivalent of the SPH equations are 
identical to the continuum alpha model, and we conjecture that they will have the 
same desirable features of the continuum model including the reduction of energy in 
the high wave number modes even when the dissipation is zero. Regardless of issues 
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concerning turbulence modelling, the SPH alpha model is a powerful extension of 
the XSPH algorithm which reduces disorder at short length scales and retains the 
constants of the motion. The SPH alpha model is simple to implement. 
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1 Introduction 



Many SPH simulations of astrophysical flows have regions where the SPH particles pass 
through each other. In a continuum simulation this streaming is unphysical because 
it implies that the velocity field is locally multi-valued. The problem can be cured by 
introducing sufficient viscosity, but this has undesirable side effects, for example excessive 
decay of shear. For this reason the version of SPH called XSPH was devised (Monaghan 
1989, Monaghan 1992). The idea is simple. Each particle is moved with an smoothed 
velocity which is constructed by averaging the velocities of its neighbours. In this way the 
counterstreaming is prevented or greatly reduced The procedure is non-dissipative and 
conserves linear and angular momentum. The smoothing has the further advantage of 
reducing local disorder. Its defect is that energy is not conserved. 

The way to conserve energy, and at the same time greatly enlarge the application of 
SPH is provided by the alpha model of turbulence (for early references see Holm, 1999). 
It is a minimalist model of turbulence which aims to take into account the finite resolution 
of any numerical simulation or experimental measurement by using an average velocity 
to move the fluid, while at the same time preserving the invariants of the flow. The latter 
is achieved by deriving the equations of motion from a Lagrangian which automatically 
provides the stress terms which are guessed in other models of turbulence (a good modern 
discussion of these models is given by Pope (2000)). The alpha equations preserve the 
circulation of the fluid (when this is an invariant of the original equations) as well as the 
additive integrals of energy and momentum. 

A feature of the continuum alpha model is that the presence of velocity gradient terms 
in the force greatly reduces the deposit of energy in the high wave number modes. In 
the SPH version the force on a particle contains non linear terms which depend on the 
velocities of neighbouring particles. We conjecture that these terms drive the energy back 
and forth from the low and high wave numbers. A related example is the Fermi-Pasta- 
Ulam problem (Fermi et al. 1955, Segre 1965) where non linear force terms lead not to 
statistical equilibrium (which would mean most energy at the high wave numbers), but 
to a flow of energy back and forth amongst a few of the modes with low wave number. A 
similar back and forth flow occurs in three dimensional turbulence (Piomelli et al. 1991, 
Jimenez 1994 , Woodward et al. 2001, Pope 2000) although this is often obscured by 
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describing just the nett flux which is from low to high wave numbers. In fact, at any 
time, different parts of a turbulent system can have energy flow either up or down the 
wave number scale. 

The advantage of the alpha model is that, if the non dissipative system already reduces 
the energy in the high wave number modes, then it is possible to introduce weak dissipation 
to complete the task of removing energy. Other turbulence models try to do this by using 
a standard form of Navier Stokes equations with a dissipation strong enough to remove 
the energy at high wave numbers. 

In this paper we derive the SPH alpha equations from a particle Lagrangian using the 
continuum constant density equations as a guide. We show that the system conserves a 
discrete form of the Kelvin circulation theorem in addition to the additive integrals of 
energy and momentum. The equations are simple to implement and the extra computa- 
tional cost is negligible. Because the turbulence model only requires minimal dissipation, 
it should be possible to reduce the normal SPH viscosity significantly in regions of the 
flow away from shocks. By combining this with the idea of varying the viscosity coefficient 
for each particle separately (Morris and Monaghan 1997) the viscosity could be reduced 
by a factor ~ 100 away from shocks. 

The plan of this paper is as follows. We derive a standard set of equations for SPH 
from a Lagrangian with a resolution length which is consistent with the density. The 
momentum equation is found to take the same form as that obtained by Springel and 
Hernquist (2001) using Lagrange multipliers to constrain the smoothing length. However, 
the thermal energy equation differs from any of their thermal energy equations. We show 
that the equations of motion lead to a discrete Kelvin circulation theorem in addition 
to the additive integrals of energy and momentum. We then define the XSPH velocity 
averaging and show that it agrees with the continuum alpha model. Guided by the 
alpha model we propose a discrete SPH Lagrangian. The Lagrangian equations of motion 
agree in the continuum limit with the equations of the alpha model. The conservation 
of circulation takes exactly the same form as in the case when the velocities are not 
averaged. Finally we discuss the dissipation and a scaling argument which suggests that 
the energy spectrum will fall off more rapidly than the Kolmogorov spectrum at short 
length scales. The extensive numerical work required to confirm the conjectures and the 
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scaling argument will be presented elsewhere. 



2 The Comparison equations 
2.1 The momentum equation 

In this section we derive a standard set of SPH equations in the absence of velocity 
averaging to provide a comparison with the alpha model equations. This standard set 
of SPH equations will be obtained for a Lagrangian allowing fully for the variation of 
resolution length h with density. The equations are similar to those of Springel and 
Hernquist (2001), and Bonet (1999, 2001) however the energy equation differs from those 
considered by Springel and Hernquist (2001). The formulation can be easily extended to 
the relativistic case using the analysis of Monaghan and Price (2001). 
The Lagrangian for compressible, non dissipative flow is (Eckart 1960) 

L = J P {^v 2 - u(p, s)) dr, (2.1) 

where u(p, s) is the thermal energy per unit mass which is a function of density p and 
entropy s. The SPH form of (2.1) with self gravity included is 

£ = E - "(ft, *) + f £ jj^jj) , (22) 

where 

if = v - (2 ' 3) 

In these equations mj is the mass, V& the velocity, r& the position, s& the entropy, and p b 
the density of SPH particle b. In practice the gravitational term is smoothed to remove 
the singularity when r b = r k . 

In the SPH formulation the density of particle a can be written 

Pa = 1 £m b W ab (h a ), (2.4) 

b 

where h a is assumed to be a specified function of p a which we denote by H{p a ) or H a . In 
many astrophysical calculations H a oc (1/ p l J d ) where the number of dimensions is d, but 
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a more general function could be used. For example, to prevent arbitarily large h when p 
becomes very small we could choose 

A 



l + Bpl /d ' 



where A and B are constants. Furthermore, while the usual practice is to estimate p a 
at a given time using the value of h a from a previous time, it would be possible and 
more consistent to calculate p a self consistently from (2.4) as suggested recently by Bonet 
(2001). Equation (2.4) is a function defining the scalar quantity p a and it can be solved 
by any standard root solving algorithm using either the value at the previous time step 
as a starting value, or by estimating the value using the rate of change of the density. 

The equations of motion follow from varying the action keeping the entropy constant. 
From Lagrange's equations for particle a 

d ' dL ^ £=0, (2.5) 



dt \ <9v a / dr n 



we find 



From (2.4) 



dv a ^ (du\ dp b „ v ^m 6 (r a -r fe ) 



^JT = £ m c V a WUK)5 cb - m a W b Wa(h), (2.7) 

where V a denotes the gradient of W taken with respect to the coordinates of particle a 
keeping h constant (we keep this convention throughout this paper), and 

6 = 1 _ H ^ mc ^). (2 . 8 ) 

where H' denotes dH b /dpb- If the density variation is smooth then £1 = 1 + 0(h 2 ). 
The first law of thermodynamics, when the entropy is constant, gives 

(du\ P . . 

where P is the pressure. The acceleration equation (2.6) with (2.7) and (2.9) then becomes 
dv n 



*- = -?"* (4 V "^ (W + ^k VM) ) ~ G ? l^F' (2 ' 10) 

The equation of motion (2.10) is exactly the same as the equation of motion due to 
Springel and Hernquist (2001) who introduce constraints on h with a Lagrange multi- 
plier. If the fluid is in a container the container can be represented by boundary forces 
(Monaghan 1994, Monaghan and Kos 1999). 
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2.2 Conservation Laws 



The symmetry of the Lagrangian leads immediately to the conservation laws. In particular 
the density (and therefore P when s is constant) and the gravitational potential are 
independent of translations and rotations of the coordinate system. Accordingly, linear 
and angular momentum are conserved. If the system is viscous then the viscous term can 
be constructed so that the linear and angular momentum are still conserved (see §2.3). 
For the constant entropy case the conservation of energy can be shown directly by noting 
that 



If the entropy remains constant during the motion the particle system is invariant to 
other transformations. Consider, for example a set of particles each with the same mass 
and entropy marking out a necklace. Imagine each particle in the loop being shifted to 
its neighbour's position (in the same sense around the necklace) and given its neighbour's 
velocity. Since the entropy is constant, nothing has changed, and the Lagrangian is there- 
fore invariant to this transformation. A discussion of the resulting invariant quantity has 
been given by Monaghan and Price (2001). We repeat the key elements of the argument 
here because we intend to show (see §6) that despite moving the fluid with a smoothed 
velocity field the same invariant holds. 

The change in L can be approximated by 



where j denotes the label of a particle on the loop. Because the system is invariant to 
the transformation we take 51 = 0. The change in position and velocity are given by 




(2.11) 




(2.12) 



5r j = r j+1 



(2.13) 



and 




(2.14) 



Using Lagrange's equations we can rewrite (2.12) in the form 




(2.15) 
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Recalling that the particle masses are assumed identical, we deduce that 

|E v ,-(r,+i-r,) = 0, (2.16) 

so that 

C = J2v r (r J+1 -r 3 ), (2.17) 

j 

is conserved to this approximation, for every loop. The conservation is only approximate 
because the change to the Lagrangian is discrete, and only approximated by the first 
order terms. However, if the particles are sufficiently close together (2.17) approximates 
the circulation theorem to arbitrary accuracy. A related argument was used by Feynman 
(1957) to establish from the invariance of the wave function that circulation should be 
quantised. 

The system is also invariant to the particles shifting around the loop in the opposite 
sense. This gives an approximation to the circulation with the opposite sign to that 
above. If these two are combined (taking account of their signs so we subtract one from 
the other) we get 

l^.toziU. (2 , 8) 

which is a better approximation to the circulation of the continuous fluid. The errors in 
the discrete approximation to the circulation theorem arise because of the higher order 
terms in the change to the Lagrangian. These errors become less as the spatial scale of 
variation of the velocity field becomes greater than the particle separation. If the velocity 
is constant then 

C = h ■ J2(v j+1 - r^) = 0, (2.19) 
as expected. If the velocity is due to rigid rotation with angular velocity fl, then 

c = £n-E r i x te+i- r i-i)- ( 2 - 2 °) 

Taking the origin of the coordinate system within the loop for convenience, the area of 
the triangle made by r,- and is \{vj x ij+i). The area is counted twice so that if the 
loop lies in a plane with normal n 

C = 2fJ-nA, (2.21) 
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where A is the area within the loop of particles. The vorticity is 2ft and (4.18) gives the 
same result as would be obtained by applying Stokes theorem to the circulation integral. 
Further examples are given by Monaghan and Price (2001). 

Salmon (1988), following Bretherton (1970) also establishes the circulation law by 
appealing to the invariance to particle interchange. However, because his analysis is 
based on the continuum equations, it is more complicated than our derivation. 

2.3 Comments on thermal energy 

If the entropy is constant we can obtain the rate of change of thermal energy from the 
rate of change of density. The rate of change of density is obtained by differentiating (2.4) 
with respect to time. We find 

dj jr = ^ J2 m ^a b -^ a W ab (h a ), (2.22) 

CLt iia ^ 

where v ab = (v — v&). The rate of change of thermal energy when the entropy is constant 
is then given by 

which disagrees with all of the thermal energy equations proposed by Springel and Hern- 
quist (2001). 

Our thermal energy equation can also be deduced by first forming the vector dot 
product of m a w a and (2.10), followed by summation over a. We get 



dE 



K 



dt 



= ~ E E W»«v a • ( V a W ab (h a ) + J^V a W ab (h b )) - (2.24) 

a b \ttaP 2 a ^bP b J dt 



where Ek is the kinetic energy and Eq is the gravitational energy. If the labels a and b 
are interchanged in that part of the double summation involving P b we can write (2.24) 

as 

— + -it = -E m «7T^E m ^ ' V a Wab{h a ) (2.25) 

dt dt Y n ap 2 ab 

Conservation of total energy then requires that the right hand side of (2.14) is —dE th /dt. 
Thus 

— - = E m "~7T = E m "7T^ E m b v ab ■ V a W ab (h a ) (2.26) 
at a at a "aP a b 
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from which (2.23) follows. 

In SPH calculations viscous dissipation is included by adding 

\^ab\ 

to the pressure terms (Monaghan 1997). The quantity o a b is a positive definite parameter, 
invariant to the interchange of a and b. It is related to the actual viscosity coefficient or 
to an artificial viscosity coefficient dependent on the speed of sound (see Monaghan 1997 
for the particular form given in (2.27), and Monaghan 1992 for an alternative IT ab ). 
Accordingly, we add the term 

-W^abVJVab, (2.28) 

1 b 

to the right hand side of (2.10) where 



VaW ah = -(V a W ab (h a ) + V a W ab (h)). (2.29) 

Angular and linear momentum are still conserved because the viscous force is along the 
line joining the centres of the SPH particles. Following the argument leading to (2.26) we 
find 

du P 1 

-IT = -^T E m bVab • V a W ab {h a ) - - ]T m b Ii ab w ab SJ a W ah . (2.30) 



We can write V a W ab in the form r ab F ab , where F ab < 0. It then follows that the 
viscous contribution to du/dt is > 0. It is also easy to show that when the fluid is viscous 
the entropy change is > 0. 

There is another approach to the effect of dissipation which is closely analogous to that 
used in the simulation of liquids at normal temperatures and pressures. In these cases 
the entropy change can usually be neglected and the dominant effect of the dissipation 
is to reduce the kinetic energy to zero. This is the case for water initially sloshing in a 
laboratory tank. In this case we write (2.24), with the dissipation term added, in the 
form 

dE* d 1 

= (E k + E G + E th ) = -- E rn b Tl ab v ab Vjy a b, (2.31) 
dt dt ^6 

where only the pressure terms are included in the rate of change of thermal energy. 

Because the right hand side of (2.31)is negative definite the SPH equations guarantee that 



10 



dissipation will cause E* to decrease as it should. For example, an oscillating polytrope 
simulated keeping the entropy of every particle constant, will evolve to a state where E* 
is a minimum. When it reaches this state Ek will be zero. 



3 Averaging the velocity field 

In the standard SPH calculations the the position of particle a is found by integrating 

£=v.. (3.1) 

However, it is possible for particles to stream through each other unless the viscosity is 
large. For that reason the XSPH version of SPH was suggested (Monaghan 1989) where 
the particles were moved with with a smoothed velocity, but the acceleration equation 
was unchanged. The smoothed velocity v a is defined by an average over the velocities of 
neighbouring particles according to 

v a = v a + e^2—(v b -v a )K(a,b), (3.2) 

b Pab 

where e is a parameter that is typically chosen to be 0.5, p ab denotes an average density, 
for example 

Pab 2 \Pa Pa) 

and K(a, b) = K(b, a) is a smoothing kernel. The kernel K(a, b) does not have to be the 
same as the kernel used in other SPH equations. In this paper we assume the kernel has 
the form 

K(a, b)= l - (W ab (h a ) + W ab (h b )) . (3.4) 

The SPH particles are moved according to 

dr. 



■ a. 



it *- < 3 ' 5 > 
The smoothed velocity has a number of useful properties. First we note that 

J2 ™ a V a = J2 m aVa, (3.6) 

a a 

which shows that the linear momentum is the same whether computed with v or v. 
Second the centre of mass moves with the v because 

m « r « = m a^a- (3.7) 
flJ a a 
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Third, the angular momentum is unaffected because 

J2 m a T a X ^ = ^ J2 m « r « X V "" ( 3 - 8 ) 

Finally we note that the motion can be exactly reversed by reversing v and v. 
The continuum limit of (3.2) is 

v a = v a + -i- / (p(r) + p a )(v(r') - v(r a ))X(r', r a )dr', (3.9) 

where 

K(r', r a ) = i(W(r - r a , /i a ) + W(r - r a , %(r')) (3.10) 
Taylor expansion about r a followed by integration gives the approximation 

a 2 1 8H do 2 

vl = vl + -V- {po>) + -(Vp • W)^^-, (3.11) 
Pa 2 dp a dh a 

where 

a 2 = € -J(x'- x a ) 2 W(r' - r a , h a )dr', (3.12) 

defines the a in the alpha model although, in our analysis, a may be different for each 
particle and, because h a can vary with time, so can a. If p is constant (3.11) reduces to 

vi = vi + a 2 (V 2 v>) a , (3.13) 

as in the alpha model. 

For typical SPH kernels a 2 oc h 2 . Accordingly, provided the scale of variation of 
the velocity field is much greater than h, the approximation (3.11) will be satisfactory. 
However, our original XSPH smoothing (3.2) applies regardless of the scale of the velocity 
field variation. In particular we expect that it can be used even when shocks occur. 

The results of this section show that the XSPH smoothing has a number of attractive 
features. It prevents SPH particles passing through each other, it retains the conservation 
of linear and angular momentum, and it agrees with the alpha model when the density is 
constant. We now show that the equations of motion can be obtained from a Lagrangian 
which is a function r and v. 
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4 Lagrange's equations for XSPH 



The analysis of Holm (1999) suggests that the Lagrangian will take the form 

L = Y > m h (^v b .y b -u{p b ,s b )), (4.1) 

where the canonical coordinates for particle b are r b with rate of change v 6 . In the interests 
of simplicity we have not included gravitational forces. These can be easily included in 
the final equations. 

The kinetic energy part of this Lagrangian can be rewritten by inverting (3.2) to give 

v a = v a -eY / —(v b -v a )W ab , (4.2) 

6 Pab 

where the errors in the inversion are of 0(h 4 ). We then get 

i^m fe v b - v 5 = ^ m 6 v 6 • v 6 + 7 E E m f nk (y b - v k ) 2 W kb (4.3) 

1 b 1 b 4 k b Pbk 

which is positive definite. 

If the sums over k are approximated by integrals, and we assume the density is con- 
stant, Taylor series expansion of the velocity difference in the integrand, followed by 
integration, shows that 

E — (v 6 - v a ) 2 W kb = (W 6 • VSj) / (x b - xfW{v b - T>)dT>, (4.4) 

k P^b J 

where v l a denotes the i th component of v Q , and there is an implied summation over i. 
Substituting this result into (4.3) gives 

lj2 m b(tf + a 2 Vvl-V.vi). (4.5) 
/ b 

In the absence of dissipation this quantity is an invariant for an incompressible fluid 
(Holm 1999) and is the alpha equivalent of the kinetic energy. Note particularly the 
terms which play a key role in the distribution of energy. The SPH equivalent is the 
sum involving the square of the velocity differences. This sum only involves neighbour 
particles, and is small if they have similar velocities. It becomes large only when the 
velocities of the neighbouring particles differ in sign or magnitude or both. This occurs 
when there is substantial change in v on the scale of h. Such a change is equivalent to 
the existence of high wave number modes. 

13 



Returning to Lagrange's equations the canonical momentum for particle a is 

-^r- = m a v a + - VV- — (v 6 - v k )(5 ba - 5 ka )W kb , (4.6) 

which reduces to 

dL ^ m a m k ^ 

-rpr- = rn a v a -e}_^— (v fc - v a )W ka , 4.7 

dv a V Pa* 

= rn a v a (4.8) 

so that the canonical momentum of the Lagrangian defined with v is the normal momen- 
tum constructed with v. This dual relation between v and v reminds us of the relationship 
between Eulerian and Lagrangian velocities noted by Holm (1999). 

The remaining term required for the equation of motion of particle a is dL/dr a . The 
details of this calculation are given in the Appendix 1. We find 

£ = " E ™» ((^ - \ - j£i + ) + (« - »)) • (4-9) 

where (a — > b) denotes terms that are identical to the previous but a and b are interchanged 
except in V , and the other quantities, where not previously defined, have the following 
meaning 

Vlb = K - Vfe) 2 , 

Ca = \ T, m A(W ak (h a ) + W ak (h k )), 
k 



and 



k Pak dh a 

The equations required for the SPH alpha model are the acceleration equation (5.9), 
and the density equation either in the summation form (4.2) or the continuity equation 

^ = 7^£™ 6 v a6 .V a W a6 , (4.10) 
at U a ^ 

together with dr a /dt = v a . The derivative following the motion is given by 

d 8 

l = ^ + v-V. (4.11) 

dt dt y ' 

These beautifully simple equations comprise the SPH alpha model. They contain all 
the desirable features of Holm's (1999) incompressible flow equations but generalised to 
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compressible flow. If the fluid is self gravitating then the gravity terms can be added as in 
§2. These new equations differ from those in §2, for example (2.10), because the velocity 
averaging leads to velocity dependent terms in the force and the particles are moved with 
the smoothed velocity. Otherwise the equations are very similar. 

In order to compare our results with that of the incompressible alpha model we assume 
h and p are constant. In this case Q = 1 and the term involving v vanishes because Vh 
is zero. The equation of motion becomes 

^-?-(S-i(|-|)) v » w + ( - 6) )' (4 - 12) 

If the sums in the SPH equations are converted to integrals the dominant terms give 
the continuum alpha model (for the details see Appendix 2). In particular, we recover 
Holm's acceleration equation (Holm 1999, equation (143)) 

dv j , e dv t 1<9/ 1_ _ a 2 dv k dv k \ , 110 , 

_ + (v . VK + v'g- = ~ w [P - -v • v - - ) . (4.13) 

We have therefore achieved our aim of constructing a Lagrangian based turbulence 
model which reduces to the alpha model in the continuum limit. In the following sections 
we will investigate some of the properties of the SPH model. 



5 Conservation laws 

In the absence of boundaries and external forces the Lagrangian is invariant to translations 
and rotations. Linear and angular momentum are therefore conserved. Because the 
Lagrangian has no explicit time dependence there is an invariant 

f)T 

E = Ei-i-i, (5-1) 

= H m « Qv • v + u(p, s)) . (5.2) 

which can be appropriately called the energy. Using (5.3) we can write 

E = ^2m b I iv b • v 6 + ^^2—{v b -v k ) 2 W kb + u{p b ,s b ) ) . (5.3) 
b \ 2 4 k Pbk ) 

If the entropy is constant then the system is also invariant to the necklace transforma- 
tion considered earlier (§2.2. In the present case the particles, assumed to have the same 
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mass and entropy, are shifted to the adjacent position on the loop and given the v of the 
new position. 

The change in L can be written 

^ = Eg ^ + |^4 (5.4) 

where j denotes the label of a particle on the loop. The change in position and velocity 
are given by 

Sr J = r i+i - r j, ( 5 - 5 ) 

and 

= Vj+i - Vj. (5.6) 
Using Lagrange's equations we can rewrite (5.4) in the form 

i(E -*))=<>. (57) 

and recalling that the particle masses are identical, we deduce 

^E v i-(ri+i-r i ) = 0. (5.8) 

which is identical to the discrete circulation theorem (2.16). The reader will note that 
this result depends on the fact that the canonical momentum of particle j is rrijVj. 

6 Is there a steady state ? 

The SPH alpha model, in the absence of dissipation, is a conservative mechanical system 
derived from a Lagrangian. If the system is at rest, for example in a uniform periodic box, 
then small disturbances would propagate as sound waves. In the linear approximation the 
waves are not coupled so that energy put into a mode stays there. If the motion becomes 
non linear the waves are coupled and any is shared amongst the modes. If equilibrium 
statistical mechanics was valid for this system each mode would have the same energy, 
and the distribution of energy between wave number k and k + dk would be proportional 
to k 2 in three dimensions. However, it is known from other cases, for example the Fermi- 
Past-Ulam study of coupled non linear oscillators (Fermi et. al 1955, Segre), that a non 

16 



linear non-dissipative system may not evolve to statistical equilibrium, but instead moves 
energy back and forth between a relatively small number of modes. 

In the present case we have non-linear coupling through both the velocity terms in the 
equation of motion, and the density variations with velocity. The results of Chen et. al 
(1999) for the continuum alpha model suggest that in wave number space the non-linear 
terms act to turn the spectrum down for large wave numbers. It is as if velocity terms in 
the equation of motion acted as a potential preventing a piling up of energy at the high 
wave numbers. This suggests that the non-dissipative model transfers energy back and 
forth between high and low wave numbers. 

The importance of this for turbulence is the following. It is known that a direct 
numerical simulation without dissipation leads to a piling up of energy in the high wave 
number modes. If dissipation is included, but it is very weak at the resolution of the 
numerical simulation, then the piling up will still take place, though eventually, if there is 
no input of energy it will decay. If there is an input of energy at the low wave numbers, 
and it is sufficiently large, then the piling up will continue because it is arriving faster 
than it can leave. The idea of sub-grid models is to prevent this piling up by including 
extra dissipation. The alpha model has the ability, without dissipation, to reduce the 
piling up of energy at high wave numbers. For this to work it is necessary for the system 
to transfer energy back and forth between high and low wave numbers. As noted earlier, 
in three dimensional turbulence there is transfer both ways (Piomelli et al. 1991, Jimenez 
1994, Pope 2000, Woodward et al. 2001) and, at any time, parts of the domain may have 
transfer up, and parts may have transfer down. The downward transfer is the residual 
between two much larger oppositely directed fluxes. The alpha model appears to mimic 
these features. How this occurs requires a numerical investigation of the non-dissipative 
system and this will be given elsewhere. 

If energy is put into the system then it is still necessary to introduce dissipation but 
it can be very much weaker than the dissipation used in the usual sub-grid models. 

7 Dissipation 

The previous analysis concerns the non-dissipative fluid. If it is correct that the non- 
dissipative SPH model transfers energy between low and high wave numbers to prevent 
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the piling up of energy at high wave numbers, then it should be possible to add a small 
amount of dissipation to remove energy at the resolution scale (~ h in the SPH calculation) 
so that, with an energy input at low wave numbers, the system reproduces the Kolmogorov 
spectrum over all but the high wave number modes. 

The dissipation term for SPH is based on Tl ab (see (2.27)). In the present case where 
the particles are moved with v we replace the velocity v in U ab with v, and substitute the 
term in the right hand side of (4.9) to obtain the dissipative equations. For convenience 
we define 



L K = 12 m a\v a ■ v , 



then (4.9) can be written 

d (dL K \ _ dL 
dt \ <9v a / <9r 



— +m a J2 m b ( Jr^VaWabiha) + (d -> 6) j = -77l a 12 m b Tl ab V a W ab . 

(7.1) 



We can deduce the thermal energy equation by requiring the total energy to be constant. 
Using the same arguments as in §2.3, but now taking the dot product of the acceleration 
equation with v a , we get the following equation for the thermal energy of particle 

d u p a 1 

-7T = ~2FT 12 m bVab ■ V a W ab (h a ) - ~12 m b U ab v ab ■ V a W ab . (7.2) 

Pa^^a b * b 

The term involving the dissipation is positive definite so that the SPH dissipation increases 
the thermal energy as it should. It is easy to deduce that the dissipation also guarantees 
that the entropy increases for each particle. 

In the same way E* , the quasi energy appropriate for systems computed assuming the 
entropy of each particle is constant (see §2.3), changes with time according to 

dE* 

—jT = 1212 m am b Tl ab v a ■ V a W ab , (7.3) 

so that, as in the case without velocity averaging, the system will evolve to a state with 
E* a minimum. 

In order to establish these results it has been necessary to use the smoothed velocity 
v instead of v in the viscosity term Ii ab . We can therefore expect that the SPH viscosity 
will be smaller than usual in those parts of the flow where the velocity is smooth. In 
SPH calculations the viscosity term n has a coefficient (also denoted by a) which can be 
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assigned a value ~ 1. If desired, each SPH particle can have its own a which can then be 
determined by the dynamics (Morris and Monaghan 1997) in such a way that its value 
is ~ 0.1 away from shocks. If the smoothed velocity is used in Ii ab) and the coefficient 
is varied as described by Morris and Monaghan (1997), then the coefficient away from 
shocks can be expected to be reduced and values of ~ 0.01 appear possible. Preliminary 
investigations of polytrope oscillations and Kelvin-Helmholtz instabilities show that this 
conjecture is true in these cases. A full investigation of these effects will be given elsewhere. 

8 Scaling and Turbulence 

The previous sections have laid out the formalism of the SPH alpha model of turbulence. If 
there is no dissipation we have argued that there is a back and forth flow of energy between 
low and high wave numbers. If there is an input of energy at low wave numbers, and the 
system has a small dissipation of the type considered above, then we conjecture that the 
SPH system will reproduce the Kolmogorov spectrum for the low wave number, energy 
containing modes. Support for this conjecture is provided by the continuum calculations 
of Chen et al. (1999). 

Some features of this process can be guessed by scaling arguments (Kolmogorov and 
Oboukov as described by Landau and Lifshitz (1993) which assume (a) that the details 
of the dissipation are irrelevant to the low wave number modes and (b) the rate of loss of 
energy through the various length scales is constant. 

We start with the energy 



and assume there is a specified rate of loss of energy from the large scale motions 8 and 
this cascades through the various length scales until it is dissipated by either by our weak 
SPH viscosity or by molecular viscosity. 

The rate of change of energy per unit mass £ determined by the first term in E for 
velocity variations with characteristic length t is 



E = m h -v 6 • v 6 + - — (V6 - v k ) 2 W kb + u(p b , s b ) 
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where the term involving a 2 can be estimated from the continuum limit of the energy 
which introduces the term a 2 Vv j ■ Vv j (see (4.5)). Because a varies with the resolution 
length the following scaling formula should be interpreted in terms of the local conditions 
of the gas. 

The characteristic velocity is 

* = < 8 - 6 ' 

If £ ^> a then this reduces to the usual Kolmogorov velocity scaling V£ = £ l / 3 S 1 ^. If £ < a 

«* = <(!)'. (8 ' 7) 

The velocity therefore reduces more rapidly with I in the smoothing region. The turnover 
time is given by 

- pr-) ■ (8 - 8) 

which is approximately constant when £ <C a. This shows that all length scales have 
essentially the same turnover time when the smoothing dominates. These effects on the 
velocity scale and the turnover time are not due to dissipation. They occur because of 
the quadratic velocity terms in the acceleration equation. 

If the energy per unit mass in the wave number range k to k + dk is E{k)dk then E(k) 
has the dimensions of (length) 3 / (time) 2 . We can relate this to the energy dissipation and 
the wave number. For large length scales we expect that E(k) will vary as S r k n . We then 
find 

E(k) = ^. (8.9) 

ks 

For the short length scales we relate E(k) to that part of S which depends on the char- 
acteristic velocity and length scales. Referring to (8.5) this is 

V= V f (8-10) 

We then find 

2 

E(k) = %. (8.11) 

The energy spectrum therefore steepens for small length scales. These approximate results 
have been confirmed for the continuum case by numerical experiments (Chen et al. 1999) 
but a similar extensive study is required to confirm these arguments for the SPH alpha 
model. 
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9 Conclusions 



The results of this paper can be summarised as follows: 

1. An averaged or smoothed velocity can be constructed for each particle such that, 
if each particle moves with this velocity, the linear and angular momentum are 
unchanged from the values with the unsmoothed velocity. This part of our analysis 
is the same as in the XSPH algorithm. 

2. The smoothing of the velocity agrees with that used in the alpha model of turbulence 
in the continuum limit with the density constant. 

3. A particle Lagrangian can be constructed using as canonical variables the smoothed 
velocity v and r. This Lagrangian has the same structure as the continuum La- 
grangian in the alpha model. 

4. The SPH Lagrangian leads to equations which conserve energy, linear and angular 
momentum and satisfy the same discrete Kelvin circulation theorem as for SPH 
with no smoothing of the velocity. 

5. In the continuum limit the new SPH equations agree with the equations of the alpha 
turbulence model when the density is constant. 

In addition we conjecture that: 

1. The non dissipative case should prevent energy piling up at the highest wave num- 
bers allowed by the resolution and it does so by transferring energy back and forth 
between low and high wave numbers. This conjecture is supported by scaling ar- 
guments, by the fact that this occurs in laboratory turbulence, and by analogies to 
dynamical systems like the Fermi-Past-Ulam problem. 

2. Weak dissipation should be adequate to reproduce the Kolmogorov spectrum for 
isotropic homogeneous turbulence. This conjecture is supported by the calculations 
for the continuum case by Chen et al. (1999). 

3. The artificial viscosity used in SPH calculations can be very much smaller than 
that currently used in smooth parts of the flow provided each particle has its own 
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viscosity coefficient. This conjecture is supported by preliminary calculations of the 
oscillations of a polytrope and the growth of Kelvin Helmholtz instabilities. 

This paper has concentrated on the dynamics of the SPH system with and without 
viscous dissipation. We have had in mind those turbulent processes which arise from me- 
chanical sources such as Kelvin Helmhotz instabilities in jets. We have not considered the 
more complex case of turbulence produced by thermal effects as in turbulent convection. 
In particular we have not considered how the smoothing of the velocity field affects the 
transport of heat. The formulation of the alpha model gives us two velocity fields and 
from these it appears possible to model the increased diffusivity expected from turbulence. 
How this should be done, however, is not clear. 

Another important issue in turbulence dynamics is mixing (Warhaft 2000). The prob- 
lem is often counter-intuitive, with many seemingly laminar flows leading to significant 
mixing (Aref et al. 1989). The natural approach to mixing is the Lagrangian formulation 
and it is possible that the SPH formulation provides an efficient practical means of deter- 
mining the extent of mixing both in turbulent and in nearly laminar flow. This is another 
area of research which is opened up by the theoretical developments in this paper. 
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A The dW/dh terms 

The contribution to the Lagrangian from the thermal energy gives the usual terms involv- 
ing the pressure. The velocity terms in the Lagrangian are 

4 b k Pbk 

where 

W bk = l -{W bk {h k ) + W bk (h b )). (A.2) 

By a change of label it is easy to show that p bk can be replaced by p b in (A.l). Then 
dL v < x . x . 2 / d f 1 \ 1 dW bk \ 

Making use of the earlier expression for the density gradient (2.7), and some straightfor- 
ward simplification, the first term in (A. 3) can be written 

"XE f-^V a WUA ) + -^VaWatihj) , (A.4) 

4 V \Pa^a Pi^h J 

where 

Ca = E m ^W 7 afc. (A.5) 

k 

The second term in (A. 3) is more complicated. We first consider the contribution from 
W bk (h b ) inW bk . We find 

7 E E fv b W 6fe (5 o6 - 4 a ) + (A.6) 

with a similar term from the W bk {h k ) in W^. Combining these terms the contribution 
not involving dW/dh is 

^ E m ^ (V a W ab (h a ) + V a W ab {h b )) . (A.7) 

4 V ^ 

The dW/dh terms give the terms involving H' in (4.9). 
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B The continuum acceleration equation 

From (4.9) the first term we consider is 

e -Y, ^ b Z^ a? VaW ah . (B.l) 
2 V Pab 

To simplify the analysis we only consider the two dimensional case and assume the density 
is constant. Without loss of generality we consider only the x component. 
The continuum limit of the summation is 



- J [v(r) - v(r')] 2 —dx'dy'. (B.2) 



Expanding v(r') about r to second order, and only retaining terms of odd parity in (x' — x) 
and even parity in (y' — y), (B.2) becomes 

\ J {-^ ~ ~ + - ->) <-> 

where v x and v y denote the x and y components of v and Ax and Ay denote (x f — x) and 
(y' ~ y) respectively. The notation v x — > v y means include a set of terms identical to the 
previous but with v x replaced by v y . After integration by parts and noting that 

J{x-x') 2 Wdx'dy' = f{y - y') 2 Wdx'dy' , (B.4) 

(B.2) becomes 

2 / dv x d 2 v x dv x d 2 v x ^dv x d 2 v x , x y ,\ (B 5) 

\ dx dx 2 dx dy 2 dx dxdy J ' 

Combining terms we can write the previous expression as 

' Q v x Q v v Q 



a 



2 KJV n2.,l 



V V + 7^ V ^ + 7T ( W ■ VyX + W • VyV ) ■ ( R6 ) 



dx dx dx 

or, on replacing v x and v y by v 1 and v 2 respectively 



a 2 ( ^! V V + A(W ■ Vv e ) . (B.7) 
\ ox ox / 



with summation on repeated indices. 

The remaining contribution from (4.9) is 

e ^ ( (a_ Cb \ 



-2 m 4^- + ^r V a W ab {h). (B.8) 
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The continuum expression for ( is 

( = j p(r') [v(r) - v(r')] 2 W(r - r')dr'. (B.9) 
Expanding the velocity difference in a Taylor series as before, we get the approximation 

C = 2a 2 W-W. (B.fO) 

Furthermore, noting p is constant, 

J2m b ^V a W ab - ^ / V r ^dr'. (B.f f) 

This integral vanishes because it has odd symmetry. The remaining term is 

J2m b ^V a W ab ^^. (B.12) 
b P P 

Multiplying (B.12) by — e/4 and combining with (12.7) gives 

«2 ( ^ V V + -^-(Vv e ■ W) . (B.13) 
\ ox 2 ax / 

which, apart from the pressure term, gives the continuum form of the SPH force/mass. 
The relevant term in Holm's equation (see (4.22)) is 

e dv e d /l_ _ a 2 dv j dv j \ _ t , 

& + &U v ' v + Ta?&?)' (R14) 

which can be written 

dv^ a 2 d f dv j dv j \ , nip > 

(u - W) ^ + T^(^^J- (Rl5) 

The expression for the smoothed velocity can be inverted to give v in the form 

v = v-a 2 V 2 v, (B.16) 
when the density is constant. Substituting this expression into (B.15) we get 

a ^ dx 2 dx ( cte £ dx e ) ) ( 

Expression (B.13) the continuum limit of the SPH term and (B.17), the continuum term 
in Holm's equation, are therefore identical. 



26 



